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Simple  models  for  crystallisation  processes 


D.  de  Cogan’  and  L.R.  Martin 

School  of  Information  Systems,  University  of  East  Anglia,  Norwich  NR4  7TJ  (UK) 


ABSTRACT 

Cellular  Automaton  (CA)  modelling,  the  repeated  application  of  simple  rules  can  result  in  very  complex  behaviour  and  is 
being  increasingly  used  to  simulate  physical  processes,  This  paper  outlines  the  technique  with  special  emphasis  on  ordered 
states  and  order/disorder  transitions. 
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1.  INTRODUCTION 

At  previous  crystal  conferences  in  Zakopane  presentations  have  been  made  about  the  benefits  of  TLM,  a numerical 
modelling  technique  which  at  one  extreme  can  be  interpreted  in  terms  of  electromagnetic  theory*-2-3.  At  the  other  extreme  it 
can  be  considered  as  the  repeated  application  of  a set  of  rules;  an  example  of  a cellular  automaton  (CA).  This  paper 
introduces  some  basic  concepts  of  general  CA  modelling  in  one  and  two  dimensions.  Even  a cursory  inspection  of  the 
subje4  reveals  a bewildering  expanse  of  interesting  avenues.  However,  as  will  be  demonstrated  here,  there  has  been  much 
recent  work  on  classification  and  this  helps  to  reduce  the  apparent  problems  associated  with  disparity  of  formalism  and 
dialect.  In  addition  to  introducing  the  subject,  this  paper  presents  details  of  work  we  have  done  in  this  area  and  discusses 
some  of  the  wider  work  on  crystallinity  and  phase  transitions.  We  will  start  with  a presentation  of  some 
interesting/amusing,  but  nevertheless  useful  two-dimensional  CAs.  We  will  then  revert  to  one-dimension  and  demonstrate 
some  of  the  CA  aspects  of  one-dimensional  TLM  algorithms  for  diffusion.  We  follow  this  with  a discussion  on  Wolfram’s 
classification  system  and  consider  several  specific  cases.  Finally,  we  will  return  to  two-dimensional  CAs  and  discuss  recent 
developments  in  terms  of  the  classification  of  CAs  which  distinguishes  between  those  that  lead  to  crystalline  states,  non- 
crystalline states  and  chaotic  states. 


2.  SOME  TWO-DIMENSIONAL  CELLULAR  AUTOMATA 

2.1.  Game  of  Life 

'Life'  was  invented  by  J.H.  Conway  and  was  popularised  by  Martin  Gardner  in  the  Scientific  American5-6.  It  involves  a 
rectangular  mesh  where  grid-points  are  either  alive  (state  T)  or  dead  (state  ’O').  The  population  of  any  point  at  the  next 
generation  is  determined  only  by  the  state  of  the  8 nearest  neighbours  at  the  present  generation  (see  figure  1). 

There  are  three  possible  transitions  in  ’Life''. 

If  kP(x,y)  = 0 then  k+iP(x,y)  = 1 if  three  of  the  neighbours  are  currently  alive. 

If  kP(x,y)  = 1 then  k+iP(x,y)  = 1 if  two  or  three  neighbours  are  currently  alive. 

If  kP(x,y)  = 1 then  k+iP(x,y)  = 0 if  four  or  more  neighbours  are  alive . 

If  kP(x,y)  = 1 then  k+iP(x,y)  = 0 if  one  or  no  neighbours  are  alive. 


Birth 

Survival 

Death  (by  crowding) 
Death  (by  isolation) 
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Life  is  totally  deterministic  but  we  generally  have  to  run  a simulation  to  see  what  any  but  the  most  trivial  starting 
configurations  will  yield  as  the  simulation  progresses. 


§01 

. . . . 

Figure  1.  The  eight  nearest  neighbours  in  Conway's  'Game  of  Life' 


2.2.  Forest  fires 

The  'Forest  Fire'  CA  uses  a square  lattice  but  differs  from  the  last  case  in  that  it  has  a stochastic  element.  Grid-points  have 
three  possible  states:  Red  denotes  a burning  tree,  White  (blank)  is  either  a dead  tree  or  an  empty  site  and  Green  is  a living 
tree.  The  rules  are  as  follows: 

If  kP(x,y)  = Red  then  k+1P(x,y)  = White  (the  tree  has  burnt). 

If  kP(x,y)  = Green  then  k+iP(x,y)  = Red,  if  one  of  the  eight  neighbours  is  currently  Red 
If  kP(x,y)  = White  then  k+|P(x,y)  — » Green  with  probability  p (growth) 

If  kP(x,y)  = Green  then  k+iP(x,y)  -»  Red  with  probability/(fire  initiated  by  lightning) 

Values  of p = 0.3  and /=  6x1  O'5  give  very  convincing  results4. 

2.3.  Others 

Examples  of  CAs  with  different  levels  of  random  element  that  are  discussed  by  Chopard  and  Droz4  include  percolation, 
Ising  spin  modelling  and  road  traffic  simulation. 


3.  TLM  ALGORITHM  FOR  ONE-DIMENSIONAL  DIFFUSION  AS  A CA  PROCESS 

The  basics  of  Transmission  Line  Matrix  (TLM)  for  diffusion  can  be  found  in  a variety  of  references  which  are  summarised 
by  de  Cogan7.  It  has  two  basic  steps,  Scatter  and  Connect.  The  Scatter  process  has  two  rules  Reflect  and  Transmit.  Reflect 
multiplies  any  data  impulse  which  is  incident  on  a spatial  node  by  the  reflection  coefficient,  p and  reverses  its  direction  of 
motion.  Transmit  multiplies  any  incoming  pulse  by  the  transmission  coefficient  x = (1-p)  and  retains  the  direction  of 
motion.  These  determine  the  magnitude  and  destination  of  the  pulses  for  the  Connect  process  which  determines  the  total 
population  distribution  at  the  next  generation.  The  resulting  scatter  algebra  can  be  summarised  for  a unit  pulse  incident  at 
position  (x)  from  the  left: 

Reflect  p[(x)]  p[(x  -1)] 

Transmit  x[(x)]  -»  T[(x+1)] 

Thus,  pp[(x)]  p [(x  -1+1)],  i.e.  p:[(x)] 

The  process  is  not  commutative:  pt[(x)]  (working  outwards)  -»  pi[(x+l-l)],  i.e.  px[(x)],  while  Xp[(x)]  -»  pi[(x-l-l)]  or 
px[(x-2)].  This  is  similar  to  the  CA  which  can  be  used  to  model  the  behaviour  of  falling  sand  in  an  hour-glass.  When  p = x 
this  is  identical  to  Galton's  'Quincunx'  which  was  an  early  experimental  method  of  demonstrating  that  such  processes  yield  a 
Gaussian  distribution8. 


150 


Proc.  SPIE  Vol.  4412 


Figure  3 Displacement  of  successive  elements  in  a p,  x binary  sequence  of  length  2k  plotted  against  number  in  sequence 
(a)  for  (k  =2,  4,  6),(b)  for  (k  = 3,  5). 


Even  numbers  Odd  numbers 


Figure  3 Displacement  of  successive  elements  in  a p,  x binary  sequence  of  length  2k  plotted  as  position  of  (s+l)th  in 
sequence  (vertical  axis)  versus  position  of  (s)th  (horizontal  axis)  for  even  and  odd  length  sequences. 

If  we  assume  that  the  diffusion  process  described  by  TLM  is  an  arrangement  by  which  data-pulses,  starting  at  the  source 
take  all  possible  paths  then  we  can  describe  the  situation  three  time-steps  after  the  start  of  our  simulation  by  means  of  a 
three-bit  binary  sequence  of  scattering  operators  ppp[(x)]  , ppx[(x)] , pxp[(x)] , pTX[(x)] , Tflp[(x)] , xpx[(x)]  , xxp((x)]  , 

and  xxx[(x)]  . Similarly,  the  scattering  algebra  can  be  used  to  determine  the  final  position  due  to  the  sixteen  four-bit  binary 
operations  after  the  fourth  time-step.  If  we  can  take  a large  binary  sequence  and  estimate  the  final  positions  of  successive 
components  then  it  can  be  seen  from  figure  2 that  there  is  a high  degree  of  correlation.  This  can  also  be  investigated  using  a 
technique  borrowed  from  phase-space  plots  in  chaos  theory.  A plot  of  Ps+i,  the  position  of  the  (s+l)th  component  in  the 
sequence  vs  Ps,  the  position  of  the  (s)th  component  yields  some  very  interesting,  and  as  yet,  unexplained  results,  which  are 
shown  for  odd  and  even  sequences  in  figure  3. 
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4.  WOLFRAM'S  CLASSIFICATION  SYSTEM  FOR  ONE-DIMENSIONAL  CAS 

The  few  examples  cited  above  give  some  indication  of  the  complexity  of  the  subject.  Stephen  Wolfram9  made  a systematic 
study  of  one-dimensional  CAs  and  recognised  that  if  there  are  k states  involving  r nearest  neighbours  then  the  total  number 
of  distinct  rules  is: 

kk7rA 

Thus,  if  there  are  two  states  (1,0)  and  if  the  outcome  at  position  (x)  is  determined  by  two  nearest  neighbours  (x+1)  and  (x-1) 
then  there  are  2s  rules.  Wolfram  used  the  same  binary  approach  as  we  have  used  in  the  previous  section  and  these  are 
shown  for  some  rules  in  the  tables  below.  The  top  row  represents  the  current  states  of  (x-1),  (x)  and  (x+1).  The  bottom  row 
in  each  shows  the  state  at  (x)  at  the  next  time  as  a result  of  applying  the  particular  rule. 


k 

dll) 

(110) 

(101) 

(100) 

(Oil) 

(010) 

(001) 

(000) 

k+1 

0 

0 

1 

1 

0 

0 

1 

0 

If  we  take  the  sequence  001 10010  as  a binary  number  then  it  represents  decimal  50.  Wolfram  termed  this  as  'Rule  50'.  The 
exclusive  OR  rule  where  the  population  at  (x)  is  one  at  the  next  generation  only  if  (x+1)  or  (x-1)  is  one  is  given  in  the  table 
below  and  it  can  be  seen  why  this  is  termed  'Rule  90'.  Chopard  and  Droz4  discuss  how  rule  184  mimics  a surface  growth 
process. 


k 

(111) 

(110) 

(101) 

(100) 

(Oil) 

(010) 

(001) 

(000) 

k+1 

0 

1 

0 

1 

1 

0 

1 

0 

Although  there  is  the  question  of  non-decidability  (just  as  in  nature)  these  rules  can  be  divided  into  clearly  empirical  classes. 


1 . After  a finite  number  of  steps  the  system  tends  towards  a unique  homogeneous  state  (i.e.  tends  towards  a limit  point). 
Rule  40  is  an  example  of  this. 

2.  The  system  develops  periodic  regions  for  almost  all  initial  states  (i.e.  tends  towards  a limit  cycle).  Rule  56  is  an 
example  of  this. 

3.  Characterised  by  aperiodic  or  chaotic  patterns,  e.g.  Rule  18  (i.e.  strange  attractor). 

4.  Rules  in  this  class  yield  persistent  complex  structures  for  a large  class  of  initial  states,  the  behaviour  can  only  be 

determined  by  explicit  simulation  (e.g.  Rule  110). 


5.  MORPHOLOGICAL  STUDIES  OF  RULE  90 

We  were  interested  in  Rule  90  for  two  reasons.  A single  initial  excitation  in  the  middle  of  a string  of  zeros  yields  the 
beautiful  Sierpinski  gasket  with  its  self-similar  properties  (figure  4). 


Figure  4 Wolfram's  Rule  90  for  different  numbers  of  iterations  after  a single  initial  excitation  (Sierpinski  Gasket) 
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The  second  reason  for  our  interest  was  the  fact  that  an  excitation  consisting  of  a random  sequence  of  ones  and  zeros  seems 
to  'develop'  a level  of  order.  We  used  techniques  borrowed  from  image  processing  to  identify  the  relative  sizes  of  triangles. 
The  situation  is  very  simple  in  the  Sierpinski  gasket.  The  smallest  triangles  have  2 ones  along  each  side.  We  call  these  class 
2 triangles.  The  next  smallest  triangles  are  class  4 (they  have  4 ones  along  each  side),  the  next  is  class  8 and  then  class  16. 
There  are  no  triangles  with  odd  numbers  of  ones  along  each  side.  The  total  number  of  class  2 triangles  at  generation,  k 
(where  k>l)  is  3k~2.  Class  4 triangles  do  not  appear  until  k = 3.  Thereafter  they  increase  as  3k'3.  The  analysis  can  be 
extended  to  higher  classes.  Any  real  simulation  of  a Sierpinski  gasket  is  bounded  by  the  size  of  the  chosen  array  and  the 
boundary  conditions  which  are  applied  (zero  boundary,  unit  boundary,  mirror  boundary,  wrap-around  boundary)  will  affect 
the  population  of  triangles.  The  effect  of  bounding  can  be  to  have  the  concentration  of  triangles  of  a given  class  tend 
towards  a steady  state.  Knowing  the  relationship  between  equilibria  and  entropy  we  investigated  the  time  variation  of 
different  classes  of  triangles  which  evolve  from  random  initial  excitations.  The  model  spaces  were  large  and  we  examined 
only  a central  region  so  as  to  eliminate  boundary  effects  within  the  number  of  iterations  used.  In  no  case  were  we  able  to 
identify  a time-dependent  change  in  the  number  of  triangles  of  a given  class.  We  are  forced  to  assume  that  within  the  limits 
of  our  numerical  experiments  there  was  no  persistent  change  in  the  state  of  order  of  the  system. 


t = 0 

t=  10 

mVm 

IJ5 

-Jl 

iSrn 

■ ■-  1 

Figure  5 Population  versus  time  for  a rule  which  outputs  zero  except  for  g( 0,3),  g(l,l)  and  g(  1,2)  which  ouput  1. 


6.  CLASSIFICATIONS  OF  SIMPLE  TWO-DIMENSIONAL  CAS 

Suzudo10  has  extended  the  classification  techniques  to  two-dimensions  but  his  analysis  may  difficult  to  assimilate  on 
account  of  differences  in  formalism.  He  defines  a generalised  CA  function,/ which  acts  on  the  present  population  of  states 
ba  to  generate  k+ia(x,y),  the  value  of  a(x,y)  at  the  (k+l)th  generation. 


a(*> y)  = /[ , a(x  - 1, y),  ta(x  + 1, y),  ta(x, y - 1),  ta(x, y + 1),  ka(x,  y)] 


He  starts  by  insisting  that  /0,0,0, 0,0)  = 0,  ones  are  not  generated  from  empty  space.  He  then  defines  a function,  g,  where  it 
is  the  sum  of  the  neighbouring  population  that  is  important,  (in  the  case  of  Life  we  are  summing  over  the  eight  neighbours) 


Proc.  SPIE  Vol.  4412 


153 


Ma(x,y ) =g[ta(x,y),  ta(x+  l,y)+ta(x  - 1,3')+  ta(x,y  -1)  + ta(x,y  + 1)] 


The  condition  g(0, 1)  says  that  ka(x,y)  = 0 and  at  least  one  of  its  neighbours  = 1 . 


g(0,l)  = 1 is  a ’Deposition’  process  and  is  the  two-dimensional  equivalent  of  Wolfram’s  Rule  90.  This  is  an  important  rule, 
otherwise  occupied  cells  may  disappear  or  become  localised  (see  figure  5).  However,  there  are  other  rules  where  we  can 
get  the  opposite  effect,  namely  locallised  empty  cells. 


We  can  use  the  reverse  of  g(0,l)  = 1,  namely  g(l,3)  = 0 which  causes  empty  states  to  advance  into  occupied  territory.  This 
is  a ’Dissolution’  or  erosion  process.  A combination  of  g(0,l)  and  g(l,3)  can  be  used  to  prevent  localisation  by  occupied  or 
empty  cells.  Figure  6 is  an  example  of  such  a mixing  process  that  can  lead  to  ’crystallisation’. 


Suzudo  then  introduces  the  concepts  of  spatial  and  temporal  entropy,  which  we  will  designate  as  Ss  and  S,  which  are 
arranged  to  be  zero  for  ordered  states  and  unity  for  disordered  states.  Ss  can  be  plotted  against  S,  for  all  possible  rules  after 
a given  evolution  time  and  it  is  possible  to  distinguish  those  rules  which  lead  to  chaotic  evolutions.  They  have  large  values 
of  both  entropies.  Even  after  2000  time-steps  there  is  a significant  clustering  of  rules  with  small  St  i.e.  less  than  0.1.  This 
allows  us  to  distinguish  between  rules  that  ultimately  lead  to  a chaotic  (liquid)  state  and  those  that  will  lead  to  a solid  state. 
In  order  to  distinguish  between  ciystalline  and  amorphous  states  we  need  another  parameter. 


The  Suzudo  p parameter  is  defined  as 


(m  = number  of  rule  entries  which  update  the  cell  with  no  change) 


This  can  be  related  to  the  X parameter  of  Langton11  which  gives  the  percentage  of  non-zero  transitions  in  the  CA  transition 
table  and  is  closely  related  to  the  Wolfram  classes  rules.  The  process  of  crystallisation  really  only  applies  to  class  3 rules. 
Accordingly,  we  can  divide  the  X,  p space  into  six  classification  areas  as  shown  schematically  in  figure  7. 
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Figure  7 Crystalline/non-crystalline  transitions  in  fyt  space  (p.  can  be  considered  as  the  CA  equivalent  of  temperature). 


7.  CONCLUSIONS 

This  paper  has  given  some  indication  of  the  complexity  that  can  be  achieved  from  the  repeated  application  of  simple  rules. 
In  spite  of  the  fact  that  the  number  of  rules  is  limited  by  the  number  of  states  and  the  number  of  neighbourhood  sites  it  is 
possible  to  distinguish  different  behavioural  classes,  and  states  of  order/disorder.  Rules  which  lead  to  an  ordered  state  have 
much  in  common  with  the  normal  processes  of  crystallisation,  even  the  inclusion  of  persistent  'defects'.  There  has  been 
relatively  little  work  on  three-dimensional  CAs.  It  may  be  that  in  numerical  terms  the  prospect  appears  to  be  daunting,  but 
with  only  two  states  and  six  nearest  neighbours,  there  are  less  rules  than  in  Life.  The  challenge  to  create  a CA  which 
emulates  real  crystallisation  processes  remains. 
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